#! /usr/bin/env python
# -*- coding: utf-8 -*-

"""
Расчет грида плотности распределения точек
Под точками подразумевается нормали к плоскостям трещин
Вариант для стереограммы
"""

from math import sqrt, exp, pi
from plane import Plane
import sys
from optparse import OptionParser
def gauss(m, sigma):
    return  exp (-1* (m)**2/(2*sigma**2) )
 
 
def line_parse(line):
    line = line.expandtabs(1)
    line = [x for x in line.split(' ') if x != '' ]    
    DIP, DIR = 0, 0    
    
    try:
        DIR, DIP = int(line.pop(0)), int(line.pop(0))
    except:
        #sys.stdout.write('WARNING! Cannot parse string: %s\n' % line)
        return None   
    line = Plane(DIR, DIP)
    return line 


# Parsing of an options
optparser= OptionParser()
optparser.add_option('-s', '--sigma', dest='sigma',
                     help='size of a net cell',
                     default=5)
optparser.add_option('-r', '--raster', dest='raster',
                     help='resolution of the grid',
                     default=50)
optparser.add_option('-j', '--intensity', dest='intensity',
                     action='store_true',
                     default=False)
options, args = optparser.parse_args(sys.argv[1:])
lines = sys.stdin.read().split('\n')

rastam = float(options.raster) # amount of rasters
sigma = float(options.sigma) # sigma in degrees; defines the size of average winidow

# parsing of the data
data=[]
if options.intensity:
    # incoming data are planes with intensity
    for line in lines:
        vals = [float(x) for x in line.split(" ") if x !='']
        if len(vals) !=0:
            data.append([Plane(vals[0],vals[1]), vals[2]])
else:
    for line in lines:
        val = line_parse(line)
        if val != None: data.append(val)  
count=len(data)
planes=data


# prepare an empty grid of rasters with its attitudes
density_grid = []
for cos3 in xrange(0, 2*int(rastam)+1):
    row=[]
    for cos1 in xrange(0, 2*int(rastam)+1):
        row.append([0, None])
        c1, c3 = (cos1*1.0-rastam)/rastam, (cos3*1.0-rastam)/rastam        
        try:
            c2 = -1*sqrt(1 - c1**2 - c3**2)
            a = Plane(0,0)
            a.define_by_plane_cos(c1, c2, c3)            
            row[cos1][1] = a            
        except:            
            pass
    density_grid.append(row)

# calculate fuction value at a zero-point to normalize all the other values
#gauss0 = 1/gauss(0,sigma)

# size of a unit window
a = Plane(0,90)
b = Plane(0,90-2*sigma)
w = a.return_angle_between(b)

# processing calc
mval=0.0
planes_dbl=[]
if options.intensity:
    for vals in planes:
        pl=vals[0]
        intensity=vals[1]
        if pl.dip<2*sigma:
            planes_dbl.append(Plane(pl.dir+180, -1*pl.dip))
        for cos3 in xrange(0, 2*int(rastam)+1):
            for cos1 in xrange(0, 2*int(rastam)+1):
                if (pl.cos1 < w) and (pl.cos3 < w) and density_grid[cos3][cos1][1] != None:
                    # !!!! this is a very ugly string
                    plane_t = pl.normal(); plane_t.dir+=180; plane_t.dip=-1*plane_t.dip
                    angle = plane_t.return_angle_between(density_grid[cos3][cos1][1])
                    P = gauss(angle, sigma)*intensity
                    density_grid[cos3][cos1][0] += P
                    mval = max(density_grid[cos3][cos1][0],mval)
    for vals in planes_dbl:    
        pl=vals[0]
        intesity=vals[1]
        for cos3 in xrange(0, 2*int(rastam)+1):
            for cos1 in xrange(0, 2*int(rastam)+1):
                if (pl.cos1 < w) and (pl.cos3 < w) and density_grid[cos3][cos1][1] != None:
                    angle = pl.normal().return_angle_between(density_grid[cos3][cos1][1])
                    P = gauss(angle, sigma)*intensity
                    density_grid[cos3][cos1][0] += P
                    mval = max(density_grid[cos3][cos1][0],mval)
else:
    for pl in planes:
        if pl.dip<2*sigma:
            planes_dbl.append(Plane(pl.dir+180, -1*pl.dip))
        for cos3 in xrange(0, 2*int(rastam)+1):
            for cos1 in xrange(0, 2*int(rastam)+1):
                if (pl.cos1 < w) and (pl.cos3 < w) and density_grid[cos3][cos1][1] != None:
                    # !!!! this is a very ugly string
                    plane_t = pl.normal(); plane_t.dir+=180; plane_t.dip=-1*plane_t.dip
                    angle = plane_t.return_angle_between(density_grid[cos3][cos1][1])
                    P = gauss(angle, sigma)
                    density_grid[cos3][cos1][0] += P
                    mval = max(density_grid[cos3][cos1][0],mval)
    for pl in planes_dbl:    
        for cos3 in xrange(0, 2*int(rastam)+1):
            for cos1 in xrange(0, 2*int(rastam)+1):
                if (pl.cos1 < w) and (pl.cos3 < w) and density_grid[cos3][cos1][1] != None:
                    angle = pl.normal().return_angle_between(density_grid[cos3][cos1][1])
                    P = gauss(angle, sigma)
                    density_grid[cos3][cos1][0] += P
                    mval = max(density_grid[cos3][cos1][0],mval)
for row in density_grid:
    for cell in row:
        sys.stdout.write('%0.4f\t' % (cell[0] / mval))
    sys.stdout.write('\n')
